o 
o 

1—5 

(N 






-)— > 

a 

a 
o 
o 



(N 

o 



cd 



o 



X 



Scaling in the growth of geographically subdivided populations: 
invariant patterns from a continent- wide biological survey* 

Timothy H. Keitt,^ Luis A. N. Amaral,^ Sergey V. Buldyrev^ 
and H. Eugene Stanley^ 

^ Department of Ecology and Evolution 
State University of New York at Stony Brook, Stony Brook, NY 1 1794 USA 

^ Center for Polymer Studies and Department of Physics 
Boston University, MA 02215 USA 

February 1,2008 



Abstract 

We consider statistical patterns of variation in growth 
rates for over 400 species of breeding birds across North 
America surveyed from 1966 to 1998. We report two re- 
sults. First, the standard deviation of population growth 
rates decays as a power-law function of total population 
size with an exponent P = 0.36 ± 0.02. Second, the num- 
ber of subpopulations, measured as the number of survey 
locations with non-zero counts, scales to the 3/4-power of 
total number of birds counted in a given species. We show 
how these patterns may be related, and discuss a simple 
stochastic growth model for a geographically subdivided 
population that formalizes the relationship. We also ex- 
amine reasons that may explain why some species deviate 
from these scahng-laws. 

1 Introduction 

Perhaps one of the most intriguing patterns in ecology 
is Taylor's law ([Anderson et ai, 1982 ; Soberon and Lo- 
evinsohn, 1987; Routledge and Swartz, 1991| ; |Leps, 1993| ; 
Curnutt et ai, 1996| ; [Maurer, 1999| ). Taylor ( |1 96 1| ) was 
the first to notice that when the mean {S) of a population 
survey is plotted versus its variance 0^(5), either in space 
or time, the relationship is typically a power-law with a 
fractional exponent 



'-iS)^{S)^ 



(1) 



Taylor was originally interested in the slope of the power- 
law relationship as a scale-free measure of spatial con- 

*In press, Philosophical Transactions of the Royal Society of London, 
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tagion or dispersion — values greater (less) than one in- 
dicate spatial clustering (over-dispersion). Later, Taylor 
used both spatial and temporal scaling as a basis for com- 
parative studies of, in his words, "synoptic population dy- 
namics" across taxonomic groups (Taylor and Woiwod, 
1982; [Taylor, 19841). 



Taylor's synoptic approach is, in many respects, a 
precursor to the recent development of "macroecol- 
ogy," (Brown, 1995), a sub-discipline of ecology and bio- 
geography (Mac Arthur and Wilson, 1967) that seeks to 
identify broad patterns in species' abundance and distri- 
bution. Macroecology has largely focussed on static pat- 
terns, such as spatial relationships between abundance and 
environmental factors (Brown et ai, 1995) and relation- 



ships between metabolic energy use and geographic dis- 
tribution (Brown and Maurer, 1987). Thus, relatively few 



continental-scale macroecological studies (Maurer, 1994 



1999) have explored interactions between spatial distribu- 



tion and population variability through time. 

In this paper, we adopt Taylor's synoptic approach and 
analyze one of the most comprehensive macroecological 
data set s available, the North American Breeding Bird 
Survey ( Peterjohn, 1994 ). The data are estimates of local 
abundance (counts) for over 600 bird species recorded at 
2-3 thousand sites (routes) across North America for the 
years 1966 to 1998. Unlike many previous macroecologi- 
cal studies, our focus is on linking geographic distribution 
to population dynamics. We report two new scaling laws, 
closely related to Taylor's power-law, for these data; one 
relating variability of population time series to their mean, 
and another relating number of sites occupied to total pop- 
ulation size. In addition, we show how these patterns may 
be related, and discuss a simple stochastic growth model 



for a geographically subdivided population that formal- 
izes the relationship. 



2 Scaling of species growth rates 

Our goal is to understand temporal variation in abundance 
at the entire population level. We therefore compute time 
series of total counts for each species by summing over all 
routes surveyed in a given year (Fig, nk). (We have pre- 
viously analyzed thes e data at the individual route level, 
see Keitt and Stanley, 1998 ) The aggregated time series 
should be relatively robust to observer err ors inherent in 
the route-level data ( Kendall et al, 1996 ) since random 
under-counts and over-counts will cancel in the summa- 
tion. After performing a bias correction (see Appendix for 
details), the resulting time series are relatively free of sys- 
tematic trends, particularly in the early years of the survey 
when the number of routes was increasing rapidly (com- 
pare dashed to solid lines in Fig. |I]a). 

We choose as our measure of the magnitude of time- 
series fluctuations, the logarithm of the ratio of successive 
counts, 

'S{t + \) 



?(0 = log 



S{t) 



(2) 



where S{t) and 5(f + 1) are the total number of birds of 
a given species counted in years t and f + 1, respectively. 
This measure has several nice properties. First, any mul- 
tiplicative, time-independent sample bias cancels in the 
ratio. Second, the measure has a natural interpretation in 
terms of population demography since, in a closed pop- 
ulation, exp[g(f)] w 1 H- (per capita birth rate - per capita 
death rate). 

As shown in Fig. [l|b, the standard deviation o{g) of 
population growth rates is strongly related to the average 
total population size. The relationship follows a power- 
law 

o(g) ~ (5)-P (3) 

for over four orders of magnitude in (5'), the total count 
averaged across years. For these analyses, we are not in- 
terested in predicting a "dependent" variable from an "in- 
dependent" variable. Rather, we are interested in mod- 
eling the functional form of the interdependence among 
variables. We therefore use major-axis regression (Sokal 
and Rohlf, 1995) to estimate model parameters. Major- 
axis regression is based on computing the leading eigen- 
vector of the covariance matrix, and minimizes squared 
errors measured perpendicular to the trend line. We also 
restrict our analysis to non-zero time series of at least 
25 years in length and with a minimum average total 
count of no fewer than 5 birds per year. Using major- 
axis regression with bootstrap precision estimates, we find 



P = 0.36 ± 0.02. Taylor's exponent (replicated across 
species) is simply y= 2(1 - P) = 1.44±0.04. 

We also study the scaling properties of population 
growth rates by examining changes in the distribution of 
growth rates with increasing population size, a technique 
familiar to statistical physicists. First, we separate the ob- 
served growth rates into three bins according to the total 
count S{t) and then construct a histogram to estimate the 
conditional probabilities p{g\S). The resulting distribu- 
tions are roughly triangular in shape with the width de- 
pending on S (Fig. nk). (The triangular shape may result 
from summing over a large number of time series with 
different local variances, see Amaral et al., 1998). If the 



distributions are "self-similar" (i.e., exhibit scaling), then 
we should be able to identify a function / that rescales the 
distributi ons so that they "collapse" onto each other (S tan- 
ley, 1971). We plot the scaled quantities: 



o{S)p 



a{S) 



aiS) 



(4) 



(Fig. pk) and find that indeed the three curves do collapse 
onto each other (Fig. gp), suggesting that p{g\S) follows 
a universal scaling form 



PislS)' 



1 



/ 



o(5) ^ \aiS) 



(5) 



These results are interesting for a number of reasons. 
In statistical physics, the presence of non-trivial scaling 
is usually taken to mean that the dynamics are largely 
governed by simple geometric properties of the system 
and do not depend stron gly on detailed properties of the 
system subcomponents ( Wilson, 1983 ). Thus, it is re- 
markable that we should find strong evidence for scaling 
across such a taxonomically and ecologically diverse set 
of species as found in the North American Breeding Bird 
Survey. These results suggest that the dynamics of North 
American Breeding Birds are unexpectedly "simple" and 
depend primarily on common patterns of internal popula- 
tion structure across species ranges, rather than details of 
individual species life-histories. 



3 Spatial structure of subpopula- 
tions 

Another reason the growth-scaling results are interesting 
is the large variability of highly abundant species. Imag- 
ine the null model that each population is subdivided into 
n equally sized, independent subpopulations, and that the 
number of these subpopulations depends linearly on S. 
The expectation, according to the central limit theorem, 
is that the standard deviation in growth rates should decay 
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Figure 1: Scaling of ecological time series, (a) Time series of total counts for three species of birds (Northern Cardinal, 
Alder Flycatcher and Brown Creeper). Dashed lines show the uncorrected total number of birds counted in a given 
year. Solid lines are bias-corrected totals (see text for explanation), (b) Standard deviations of empirically observed 
population growth rates plotted against average annual total counts for 428 species of North American breeding birds. 
Results for the time series plotted in panel (a) are highlighted. Notice that the highlighted time series are less variable 
than other species in the data set with similar total population size. 



as the -1/2 power of S. The observed decrease in fluc- 
tuations is considerably slower (i.e., p < 1/2), such that 
highly abundant species are considerably more variable 
than expected under the null model. 

To account for the increased variability for highly abun- 
dant species, we require that the number of subpopula- 
tions does not scale in a simple, linear fashion with in- 
creasing S, but instead takes the form: 



■,1-a 



(6) 



with a ^ 0. Values of a > will be found, for example, 
when the "typical" size of the subpopulations also scales 
with total abundance, i.e., there is a positive relationship 
between regional and loc al abundance, a well docum ented 
pattern in macroecology ( Gaston and Lawton, 1988 ; Gas- 
ton, 1996). The positive correlation between local and 
regional abundance results in fewer subpopulations for a 
given total population size as each subpopulation accounts 
for more individuals. Again appealing to our observation 
that, under the central limit theorem, o{g) ^ rT^''^, and in 
combination with equations (^ and (^, it is strait forward 
to show that for roughly equal-sized subpopulations, the 
estimated exponents must obey 



P 



1-a 



(7) 



We do not have access to a precise estimate of the 
number of subpopulations for each species in the survey. 



However, we can use as a proxy the number of survey 
routes where a species had a non-zero count in a given 
year. To test the assertion in Eq. ^ we plot number of 
survey routes with non-zero counts hit) versus the (un- 
corrected) total count Sit) for all bird species recorded in 
the survey in 1997, excluding species seen at fewer than 
5 routes or with fewer than 5 total individuals counted. 
The data follow closely the power law dependence pre- 
dicted by Eq. (^ with an exponent a ~ 0.25 ±0.03, again 
using major-axis regression with bootstrap precision esti- 
mates (Fig. Ha). Remarkably, the estimate of a predicts 
a value of P = 0.38 ±0.02, very close to the estimate 
(P = 0.36 ± 0.02) obtained by measuring the standard de- 
viation in growth rates directly (Fig. Hp). Even more 
striking is the consistency of our estimate of a across 
years, despite large changes in the number and spatial 
distribution of sampling locations through time (Fig. 0b). 
These results directly imply that average local abundance 
(i) — S{t)/h{t), scales with total (regional) abundance ac- 
cording to (s) ^ S{t)"-. 

We can gain further insight into the organization of a 
species population in different routes by considering how 
the distribution of number of routes with non-zero counts 
depends on total counts. That is, we may quantify the or- 
ganization of the subpopulations through the conditional 
probability density p{h\S), which measures the probabil- 
ity to find a bird species with S total counts having non- 
zero counts in n distinct routes. Fig. ^suggests thatp(«|5') 
will have a peak that increases as a power law with S. As 
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Figure 2: Distribution of growth rates for different population size classes, (a) Probability density p{g\S) of the 
growth rate g for all bird species in the North American Breeding Bird Survey database. The distribution represents 
all annual growth rates observed in the 3 1 -year period 1966-1996. We show the data for three different bins of initial 
sizes. The solid lines are exponential fits to the empirical data close to the peak. The approximately triangular shape 
of the distribution may be the result of mixing Gaussians with different widths ( [Amaral et al., 1998 ). (b) Scaled 
probability density pscai = '^p{g\S) as a function of the scaled growth rate gscai = [g^ g]/*^ for all species and years in 
the survey. The values were rescaled using the measured values of g and a. All the data collapse upon the universal 

CUrVepscal=/(-|^scal|)- 



shown in Fig. ^ this is indeed the case. If the data exhibit 
scaling, we should be able to identify a universal scaling 
function h such that 
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We test the scaling hypothesis in Eq 
scaled variables: 



}■ 




(8) 


■ (I) 


by plotting 


the 


h 




(9) 



Fig. ^ shows that all curves collapse onto a single curve, 
which yields the scaling function h{u). 

4 Discussion 

Our analys is differs from Taylor's orig i nal studies (T ay- 
lor, 1961; [Taylor and Woiwod, 198|; fTaylor, 1984|) in 



an important way. Taylor was interested in comparative 
analysis and so calculated a separate exponent for each 
species. He did this by analyzing multiple samples, repli- 
cated across time or space, for each species. Here, we 
calculate exponents replicating across species. One ad- 
vantage of this approach is that we analyze the time series 
of total counts, summed over the entire survey. These to- 
tal counts are considerably more robust estimates of abun- 
dance than local counts taken at individual routes. 



Another advantage of analyzing scaling across species 
is that is allows us separate general patterns or "laws" 
that are invariant across taxonomic groups from general 
rules that may explain deviations from these laws. Our 
reasoning is that when the physical dimensions of a prob- 
lem, such as energy or material flows, or spatial popu- 
lation structure, predominate, we should observe scaling 
laws that do not depend strongly on the biological dif- 
ferences among species, but that species-specific differ- 
ences should appear as residual variation after the com- 
mon scaling laws are factored out. That so many species 
fall along a single scaling relationship describing variabil- 
ity as a function of population size (Fig. ^1-0) suggests that 
there may be universal features to the way in which North 
American breeding bird populations are subdivided spa- 
tially. We find exactly these features in the invariant, 3/4- 
power scaling of number of occupied survey routes versus 
total population size (Fig. g|-Q). 

However, not all of the variability in the data is ac- 
counted for by these scaling laws. For example, species 
with average total counts of approximately 250 individ- 
uals exhibit more than two orders of magnitude range 
in growth rate standard deviation (Fig. |l]b). We believe 
that this residual variation does reflect important aspects 
of the ecology of individual species. There is a strong 
correlation between the residuals in Fig. [l|b and the area 
of the corresponding species ranges, measure in terms of 
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Figure 3: Statistical analysis of the number of routes populated by a given bird species, (a) Double logarithmic plot 
of the number of routes with non-zero counts n{t) versus total number of birds counted 5(f) for each species observed 
in 1997. The bias correction applied to the time series in Fig. |l| is unnecessary in this case as all data come from a 
single year. The data for all species follow closely a straight line the log-log plot suggesting a power law dependence. 
From the slope of the line, we estimate 1 — a = 0.75 ±0.03. (b) We perform a similar analysis for all 31 years in 
the database and plot the exponent estimates for each of the years. Our results show that the power law dependence 
remains remarkably stable during the 31 survey year, clustering around a = 0.25. Error bars are bootstrap 95% 
confidence intervals. 



the average number of non-zero routes (T. Keitt, unpub- 
lished results). A likely explanation for this pattern is 
that fluctuations in the abundances of broadly distributed 
species will tend to average out spatially because differ- 
ent regions are influenced by geographically distinct cli- 
mate regimes. Thus, it appears that species whose life 
histories tend to produce strongly aggregated distributions 
(i.e., species that are locally common, but regionally rare) 
are the ones that fluctuate the most relative to their total 
abundance. Species that have broad spatial distributions 
(i.e., locally rare, but regionally common) are therefore 
expected to fluctuate less than similar species with more 
restricted geographic ranges. 

Ranking species in terms of their residual deviation 
from the growth-scale law (Fig. |lp) supports our hypoth- 
esis that locally common, but regionally rare species fluc- 
tuate more than expected, and vise versa. Large posi- 
tive residuals correspond to species with restricted geo- 
graphic ranges, such as Golden-cheeked Warbler (Den- 
droica chrysoparia; 2.5 times more variable), species that 
are habitat specialists and nest in large colonies, such 
as Tricolored Blackbird {Agelaius tricolor, 13.7 times 
more variable), species that breed in large groups called 
"leks", such as Greater Prairie Chicken {Tympanuchus 
cupido; 3 times more variable), and species that show 
strong local migration patterns in response to changes 
in resource availability, such as White-winged Cross- 



bill (Loxia leucoptera, 3.5 times more variable) and Red 
Crossbill (Loxia pytyopsittacus; 3.7 times more variable). 
Species that show low variability in relation to the scaling- 
law are typically solitary, territorial breeders such as 
Yellow-throated Warbler (Dendroica dominica; 2.8 times 
less variable), Prairie Falcon (Falco mexicanus; 2.6 times 
less variable). Swamp Sparrow {Melospiza georgiana; 2.5 
times less variable), Kentucky Warbler iOporornis for- 
mosus; 2.4 times less variable), and Chuck-will's-widow 
(Caprimulgus carolinensis; 2.3 times less variable). The 
important point is that had we started from a purely aute- 
cological standpoint and ignored the important physical 
dimensions of the problem (e.g., structure of geographic 
ranges), we could easily have missed key pattern in terms 
of deviations from general scaling laws. 

We should however mention several caveats. We do 
not as of yet know whether our results can be general- 
ized to include other, non-avain taxonomic groups, or to 
other continents and climate regimes. Also, despite our 
use of highly aggregated, and therefore more robust time 
series, we suspect that there remain sources of variation 
in our analysis unrelated to actual population fluctuations. 
One vexing problem is repeated local migration between 
sampled and unsampled locations (we call this "slosh- 
ing"). Even if there is no variation in the true abundance 
across years, sloshing will lead to a given individual being 
counted in some years and not others, leading to measure- 
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Figure 4: Statistical analysis of the number of routes populated by a given bird species, (a) Conditional probability 
density function p{n\S) of finding « non-zero count routes for a bird species with S total counts. To improve the 
statistics, we partition the bird species into four groups according to size, (b) To illustrate the scaling relation (g), we 
plot the scaled probability density 5'^" p(n/5''^"|5') versus the scaled number of non-zero routes ii/S^^"', combining 
data from all years. In agreement with (M), we find that the scaled data fall onto a single curve. 



ment errors in the time series. We suspect this effect is 
not a significant component of variation in most of our 
time series, but may be substantial in a few cases. (Slosh- 
ing may contribute to the extreme variability of Tricolored 
Blackbird, for example.) Additional data, such as mark- 
recapture, may be need to resolve this issue. 

There are two additional mechanisms related to our 
model for geographically subdivided populations that we 
have not discussed. We have shown how a non-linear 
dependency of the number of subpopulations versus to- 
tal population size may explain the observed deviation 
from 1/2-power scaling of population fluctuations. Our 
basic hypothesis depends on the average local abundance 
scaling with total abundance in independently fluctuat- 
ing subpopulations of roughly equal size. However, 
there are other patterns that may influence the "effective" 
number of independently fluctuating subpopulations and 
thus partially amount for the observed exponent in the 
growth-scal ing law. First, large spatial variation in local 
abundance ( Brown et al, 1995 ) could cause wide-spread 
species to fluctuate with greater magnitude than if all sub- 
populations have the same local abundance, since most of 
the variation would be driven by a few, high abundance 
sites. Second, strong spatial autocorrelation in population 
growth increments or "spatial synchrony" among fluctuat- 
ing su bpopulations (|Grenfell et al, 1998 ; Bj0rnstad et al.. 



1999; [Kendall et al, 2000| ; [Lundberg et al, 2000D may 
also cause a reduction in the effective number of indepen- 
dent subpopulations, and thus account for the increased 
magnitude of fluctuation in broadly distributed species. 



Temporal autocorrelation may act similarly to increase 
or decrease variability relative to our model. The con- 
sequences of these mechanisms need further exploration. 

A surprising result of our analysis is the, to our knowl- 
edge, previously unreported 3/4-power scaling of spatial 
distribution as a function of total population size (Fig. g). 
This result is closely related to, but not the same as, the 
"Distribution-Abundance" curve of Hanski and Gyllen- 
berg (1997) that describes the fraction of regional habitats 
occupied as a function of average local abundance. We 
do not as of yet have an explanation for why the exponent 
should take this particular value, nor why it is so consis- 
tent through years. Recently, there has been considerable 
interest in explanations for the apparent 3/4-power scal- 
ing law rela ting body mass to metabolic output ( Enquist 
et al, 1998; [Westefg/., l"999| ; podds ef a/., 200l| ; Niklas 
and Enquist, 2001). One explanation posited to explain 
3/4-power scaling is optimal structuring of a fractal trans- 
port network, such as the vascular system of plants and 
animals ( West et al, 1999[ ). This suggests an interesting 
hypothesis to explain 3/4-power scaling in our analysis: if 
the geographic ranges of species are subdivided in accord- 
ing to a particular fractal pattern, perhaps because of the 
fractal nature of the physical environment (e.g., Rinaldo 
et al, 1995), then it might lead to our observed scaling 
laws. Testing this hypothesis will require additional study. 

It is interesting to note that our results are in striking 
qualitative agreement with similar studies from a broad 
range of social systems, ranging from growth of compa- 
nies in the U.S. economy to the GDP of countries (Stan- 



ley et ai, 1996; |Lee et al, 199^; [Plerou et al, 1999|), RefcrenCCS 



suggesting that our simple model of growth may apply 
quite broadly ( Amaral et al, 1998). Our observation that 
more "specialized" birds (in terms of smaller number of 
subpopulations) fluctuate more in total number than those 
that average fluctuations over many subpopulations may 
have an interesting parallel in social organizations: those 
that specialize on a few economic activities, e.g., coun- 
tries with a single export product, may fluctuate consider- 
ably more than similarly sized organization with diverse 
economic activities, e.g., countries that produce a range 
of products. Putting all of one's eggs in a single basket, as 
the saying goes, some times leads to catastrophes, and, it 
appears, greater variability as well. 
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Appendix: Corrections applied to 
time series 

Let Sfiiv be the number of birds of species u counted at 
route V in year t. The raw total counts R,u — Y.y' Stuv con- 
tain information about the abundance of species u in year 
t as well as information about the number M and distribu- 
tion of routes surveyed in year t. The goal is to remove the 
bias in the counts Rm introduced by variation in the num- 
ber and distribution of survey routes through time. We 
do this by replacing each count St„y for a given species at 
a given route with the time average /j„,, = T^7 ' Y.t ^tuv for 
that route and species, where Ty is the number of years 
that route v was surveyed. We then construct new, sur- 
rogate time series Mm = Y]^,' jAuv whose variation only re- 
flects changes in the number and distribution of survey 
routes through time (because the same /j,„. is used in each 
year), and not any real change abundance. We can then 
generate a bias corrected time series by subtracting these 
new time series from the raw totals: 



Rru-Mr, 



-M„ 



(10) 



where M„ is the time average of M(„ for species u. The 
advantage of this approach is that survey routes added or 
removed outside a species range will not influence the cor- 
rected total, because these routes will have /jj„, = 0. 
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